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LIST OF SYMBOLS 


En glish Symbols 

A pipe cross-sectional area 

a maximum eigenvalue 

a R Rosseland mean opacity 

c speed of sound 

C p specific heat at constant pressure 
C v specific heat at constant volume 

D diameter 

D h hydraulic diameter 

e total energy per unit volume 
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h enthalpy or heat transfer coefficient 
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p pressure 
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q" conductive heat flux 

q" radiative heat flux 
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q" wall heat flux 

q t total heat flux 

q'" volumetric heat generation rate 
Q total heat flux 

Q uniform heat generation 

r radial position measured from the centerline 

R pipe radius 

gas constant 

Re Reynolds number, uL / v 

t time 

T temperature 

T s surrounding temperature 

u axial velocity, z-component 

v radial velocity, r-component 

z axial position measured from the pipe entrance 

Greek Symbols 
a void parameter 

(5 clustering parameter 

8 Delta number kAT/(q D h a) 

8 b thickness of laminar sublayer 

8 V velocity residual 

8 T temperature residual 

e internal energy 

y ratio of specific heats, C p j C v 

r) positioning parameter 
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k Boltzmann constant 

ji dynamic viscosity 

v kinematic viscosity, p/p 

p density 

a Stefan-Boltzmann constant 

a ph photon collision cross section 
t shear stress 

|© | magnitude of vorticity 

Subscripts 

1 inside wall boundary mesh point 

2 inside fluid boundary mesh point 

b bulk conditions 

c convective components 

e eddy components 

i axial mesh point locations 

j radial mesh point locations 

m molecular component 

average 

r radial locations 

radiation components 
t turbulent component 

w wall conditions 
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CHAPTER 1 
INTRODUCTION 

1.1 Motivation for the Heat Transfer Modeling 

Numerical heat transfer is a key issue in nuclear reactor analysis. For a high 
temperature energy conversion system, energy from the nuclear fission process is 
transferred to the coolant system by conduction, convection and radiation. This study is 
especially focused on numerical analysis to solve the heat transfer problems in nuclear 
reactor core design, which include the Heterogeneous Gaseous Core Reactor [HGCR], 
High Temperature Gas Cooled Reactor [HTGR] and the XNR2000 nuclear rocket. Table 
1-1 presents the applications of numerical modeling. 


Table 1-1 Applications of numerical modeling. 


Reactor 

Fluid Flow 

Heat Transfer 

Thermal Boundary 

HGCR 

UF 4 gas 

Convection & Radiation 

Gas deposits heat to wall 

HTGR 

Helium gas 

Convection 

Gas removes heat from wall 

XNR2000 

Hydrogen gas 

Convection 

Gas removes heat from wall 


A Heterogeneous Gas Core Reactor uses a gaseous fissile material as fuel for power 
generation^ 131 This allows for power generation at temperatures much higher than the 
melting point of solid fuel nuclear reactors. Power generation and power conversion at 
very high temperatures can potentially reduce the system mass and improve the specific 
impulse performance of a nuclear thermal rocket. One of the most challenging issues 
related to design and operation of gas core reactors is the containment of the high 
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temperature fissioning plasma. The wall cooling is the most important issue in design of 
gaseous core nuclear reactors. The heat transfer process involved in an ultrahigh 
temperature gas core reactor systems is characterized by the convective flow of a radiating 
gas. Uranium compound gases at pressures in the range of 10 to 40 atm are optically thick. 
At temperatures close to 3500°K which is the typical exit temperature of the reactor core in 
a more recent design, the radiative heat transfer rate in these opaque gases is higher and 
comparable with the convective heat transfer rate. 1371 Therefore, the heat transfer analysis of 
a fissioning gas must include both convective and radiative transfer. The flow and radiation 
transport equations must be solved simultaneously in order to determine the temperature 
distribution and heat transfer rate. Figure 1-1 shows the high temperature vapor reactor 
with advanced energy conversion system. 

The HTGR primary system is composed of several loops, each housed within a large 
cylinder of prestressed concrete/ 441 The flow is directed downward through the core by a 
circulator mounted above the steam generator in the cold leg. The reactor vessel and 
steam generator are connected by a short, horizontal cross duct. The coolant from the core 
exit plenum is directed laterally through the interior of the cross duct into the inlet of the 
steam generator. Coolant from the steam generator and circulator is directed laterally 
through the outer annulus of the cross duct into the core inlet plenum. The HTGR 
primary coolant flow path is illustrated in Figure 1-2. 

XNR2000 is an expander cycle nuclear rocket engine powered by a fast-spectrum 
cermet-fueled nuclear reactor that heats hydrogen to a maximum propellant 
temperature/ 41 The reactor is comprised of an outer annulus core of Mo-U0 2 prismatic 
fuel elements and a cylindrical inner core of W-U0 2 prismatic fuel elements. The 
baseline XNR2000 reactor core consists of a total of 151 prismatic fuel elements of 55 
cm in active length. The core is arranged such that an inner-outer core configuration is 
obtained with 61 inner core fuel elements and 90 outer core fuel elements. 1141 The 
purpose of inner-outer core configuration is to provide a folded flow path for the 
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hydrogen propellant flowing upward through the outer core and downward through the 
inner core. The XNR2000 coolant flow path is shown in Figure 1-3. 

1.2 Computational Method for CFD and Heat Transfer 

The development of heat transfer modeling has been a major area of research for 
several decades in the nuclear thermal hydraulics field. Most people attribute the first 
definitive Computational Fluid Dynamics (CFD) work to Richardson in 1910, who 
introduced point iterative schemes for numerically solving Laplace’s equation and 
bihamonic equation. He clearly defined the difference between problems which must be 
solved by a relaxation scheme. In 1918, Liebmann presented an improved version of 
Richadson’s method. Liebmann’ s methods used values of the dependent variable both at 
the new and old iteration level in each sweep through the computational grid. The 
beginning of modem numerical analysis is attributed to a famous paper by Courant, 
Fredrichs and lewy, [1] the CFL, frequently seen in the literature. In that paper, uniqueness 
and existence questions were addressed for the numerical solutions of partial differential 
equations. It is original source for the CFL stability requirement for the numerical 
solution of hyperbolic partial differential equations. However, limits to the above method 
still existed for steady-state and low temperature conditions at that time. 

In 1940, Southwell introduced a relaxation scheme which was used in solving fluid 
dynamic problems where an improved relaxation scheme was required. During the 
decades of 1940s and 1950s, Southwell’s methods were generally the first numerical 
techniques introduced to engineers. Allen applied Southwell’s scheme to solve the 
incompressible, viscous flow over a cylinder. They used the empirical or semi-empirical 
information, wall function for the heat transfer studies. 1461 However, this approach is too 
sensitive to the near wall grids and inaccurate at the flow separation region. Therefore, a 
highly-stretched fine grid near wall boundaries is required to solve unsteady and turbulent 
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flow problems. Unfortunately, the fine grid yields computational difficulties in terms of 
the stability limitation and the computation time. To remove the time step limitation, the 
fully-implicit numerical schemes were developed in the mid- 1 970’ s by Briley and 
McDonald. 171 

Professor John von Neumann developed his method for evaluating the stability of 
numerical methods for solving time-marching problems. O’Brien, Hyman and Kaplan 
later presented a detailed description of von Neumann method. This paper is significant 
because it presents a practical way of evaluating stability. At same time, progress was 
being made on the development of methods for both elliptic and parabolic problems. 
Peaceman and Rachford developed a new family of implicit method for parabolic and 
elliptic equations in which sweep directions were alternated and the allowed step size was 
unrestricted. 1201 But those methods cannot be used to handle the discontinuity problems, 
such as shock capture problem. It was difficult to solve transonic and supersonic fluid 
flow problems. 

Early efforts at solving flows with shock waves were part of Lax’s approach. 191 Lax 
and Wendroff introduced a method for computing flows with shocks which was a second- 
order scheme that avoided the excessive smearing of earlier approaches. MacCormack 
has devised an implicit scheme that requires only the inversion of block bidiagonal 
systems rather than block tridiagonal systems, thus yielding savings in computer time and 
storage requirements. During the past 20 years, the explicit methods were developed to 
solve the compressible Navier-Stokes equations which include the Hopscotch method, 
DuFort-Frankel method, Brailovskaya method, Allen-Cheng method, Lax-Wendroff 
method and the MacCormack method. 1171 All of above methods, except the MacCormack 
scheme, are first-order accurate so that they cannot be used to accurately compute the 
time evolution of a flow field. 1231 In addition, all of the methods have a stability 
restriction which limits the maximum time step. The allowable time step is given by the 
CFL condition, which for 2-D problem becomes 
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where a is the maximum eigenvalue. For the other schemes, analytical stability conditions 
cannot be obtained 120 ' and a numerical investigation is presented in Figure 1-4. In this 
graph, the schemes are stable in the region below the corresponding curve. Figure 1-4 
shows that the MacCormack method presents the best stability. 

The MacCormack is the most popular two-step Lax-Wendroff method for solving 
problems with shock-capturing schemes.* 8 ' This method is designed to solve time- 
dependent equations such as the complete Navier-Stokes equations without any artificial 
dissipation term or limiters. It is based on the second-order accurate explicit predictor- 
corrector method but adds an implicit procedure in the predictor-corrector sequence for 
points at which the local CFL number exceeds the stability limit. The method has been 
applied to two-dimensional internal supersonic flows, two-dimensional external flows, 
external axisymmetric flows and three-dimensional flows over a biconic body with a 
compression flap. 128 ' This scheme was applied to either the complete or thin layer forms 
of unsteady Navier-Stokes equations. 

MacCormack’ s method is one of the most efficient of the second-order schemes from 
point of view of operation count. This approach is well adopted for time-dependent 
problems; it should be inserted in a multi-grid framework. The investigations on multi- 
grid Lax-Wendroff-type schemes can be found in reference [17]. 

The important conclusion is that among all second-order viscous Lax-Wendroff 
schemes, MacCormack scheme, with flux splitted operators, presents the best stability, 
consistency and efficiency.* 25 ' 



1.3 Objective and Overview 


The main objective of this study is to develop a computational fluid dynamics and 
heat transfer model for convective, conductive and radiative heat transfer in high power 
density gas cooled and gaseous core nuclear reactors and the XNR2000 nuclear rocket 
core. To achieve this goal axisymmetric, thin-layer Navier-Stokes equations associated 
with 1-D integral approximation, Rosseland’s diffusion approximation and algebraic 2- 
layer eddy viscosity turbulence model are used to simulate heat transfer in nuclear reactor 
cores. An implicit-explicit, finite volume, MacCormack method in conjunction with 
Gauss-Seidel line iteration procedure is utilized to solve the governing equations. An 
enthalpy-rebalancing scheme is implemented to allow the convergence solutions to be 
obtained with the applications of a wall heat flux. A two-dimensional method based on 
finite element technique is used to investigate the geometric behavior of a nuclear reactor 
fuel elements. 

Chapter 2 describes MacCormack implicit-explicit model, Lomax and Baldwin’s 
two-layer algebraic turbulence model, Rosseland’s diffusion model, wall heat transfer 
model, enthalpy rebalancing scheme and the thermal conduction model. Chapter 3 
provides the thermophysical properties of dissociated hydrogen, helium and UF 4 gases. 
Chapter 4 presents the assessment and validation of the numerical model with separation 
angle and drag force calculations around a sphere suspended in a cylinder, and heat 
transfer correlation comparison. Chapter 5 presents the applications of numerical models 
in gaseous core reactors, gas-cooled reactors and the XNR2000 nuclear rocket cores. 
Finally, the conclusions of this thesis are described in Chapter 6. 
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Figure 1-2 HTGR primary coolant flow path (From Reference [44]). 
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Figure 1-3 XNR2000 coolant flow path (From Reference [14]). 
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1 - Richtmyer scheme 

2- MacCormack scheme 

3- Lerat-Peyret scheme 

4- The CFD condition <j=l 

5- Von Neumann scheme 

Figure 1-4 Comparison of stability limits for the second-order schemes (From Reference 

[ 20 ]). 



CHAPTER 2 

NUMERICAL MODELING 


2.1 Governing Equations 

Combined convective and radiative heat transfer of simultaneously developing 
turbulent flow in a smooth circular tube with nongray gases is numerically modeled. 
Based on the physical considerations, the two-dimensional axisymmetric thin-layer 
compressible Navier-Stokes equations in strong conservation law form with Rosseland’s 
approximation radiative heat transfer model and the algebraic two-layer eddy viscosity 
turbulence model are used to simulate the fluid flow and heat transfer in the gas core 
reactor. In order to derive the governing equations, the following assumption are 
made: 1 101 

1 . The axial viscous dissipation and thermal radiation are negligible compared with 
the radial viscous dissipation and thermal radiation. 

2. The axial heat conduction is negligible. 

3. The gas absorbs and emits radiation but does not scatter. 

4. The wall is a gray diffuse emitter and reflector. 

5. There is no dissociation or phase change for UF 4 gas in the high temperature 
range. 

Based on these assumptions, the time-dependent, mass-averaged and compressible 

T 28 ] 

Navier-Stokes equations in strong conservative and axisymmetric form are given as: 1 
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p uv 
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and the viscous and thermal source terms are 


G., = 


0 

du 

^Vr 

4 *^-- 1 
3 dr 3 * T r 
du 4 dv " » 

^- + -^-- 9 , -q. 


and 


(2-1) 


(2-2) 


(2-3) 


0 

_ 0 

H - 0 

• 

Yq. 

In the above formulas/ 431 u and v are the velocity components in z and r directions, 
respectively, p is the density, P is the pressure, p T is the total viscosity, q" c is the 
conductive heat flux, and e is the total energy per unit volume which is related to the 
internal energy e and kinetic energy. 


(2-4) 
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e = p 


+ ^(u 2 + v 2 ) 


(2-5) 


The total viscosity ji is given as the sum of its molecular component p m and 
eddy component p e : 


P = P* + P e 


( 2 - 6 ) 


The eddy viscosity is obtained from the turbulent eddy viscosity model proposed by 
Baldwin and Lomax. [5] The heat flux term includes the conductive heat flux term, q c , 
and radiative heat flux term, q r . Based on Fourier’s law, the conductive heat flux can 
be written as: 



where 


(2-7) 


k=C, 


\P PJ 


(2-8) 


Typical values of the molecular and turbulent Prandtl number at standard condition 
are 0.9 and 0.72, respectively. The radiative heat flux term is presented in a later section. 

The complete set of Navier-Stokes equations (2-1) includes one continuity 
equation, two momentum equations, and one energy equation. In order to get six 
unknowns, p (density), u (axial velocity component), v(radial velocity component), e 
(total energy per unit volume), P (pressure) and T (temperature), two more supplement 
equations which are given by perfect gas law are required. Under perfect gas law, there 
are the following relations: 
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P= p RT s =C V T h = C P T y = 


C, 


(2-9) 


when using above relations, two additional equations are obtained as follows: 


P = ( y-i 





( 2 - 10 ) 


T = 


(y -0 




e l/ 2 
— 

L p 2 v 


(2-11) 


Hence, combining above equations (2-1), (2-10) and (2-11), a set of governing 
equations is formed. 

The Navier-Stokes equations are the coupled nonlinear mixed hyperbolic-parabolic 
system of partial differential equations. At high Reynolds number condition, these 
equations become stiff and are difficult to solve. Reynolds number is a measure of the 
ratio of the inertial to the viscous forces of a fluid. The viscous terms which cause the 
system to be parabolic are of the order of the reciprocal of the Reynolds number in 
magnitude. At high Reynolds number the system is almost everywhere hyperbolic, the 
viscous terms are negligible except in the boundary layers. Within these layers, the 
viscous terms are significant and control the important phenomenon of boundary layer 
separation. The disparity in magnitude at high Reynolds number between the inertial and 
viscous terms causes the system to be mathematically stiff. The only way to solve the 
complete Navier-Stokes equations is numerical procedure which will be discussed in the 
later sections. 


2.2 Numerical Procedure 


A hybrid implicit-explicit numerical method based on the finite volume 
discretization approach is used for computer implementation. The finite volume method 
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uses the integral forms of conservation equations on the finite cells. The equations are 
approximated by summing the fluxes of mass, momentum and energy from neighboring 
cells into each cell, thereby partial differential equations are changed to the discrete 
equations. 

To derive an implicit form, equation (2-1) is differentiated with respect to time , t, 


(8U_ 

'l at J 

dt 


dA, 




'dlY 

\dt > 


8B, 


dz 




diT 

dt > 


dB, 


du_ 

V dt J dH 


dr 


dr 


8t 


( 2 - 12 ) 


where A i? Bj and B v are the Jacobians of Fj, Gj and G v with respect to U, they are. 


A .- d A 

‘ dU 


(2-13) 


dG, 

B = — L 
' dU 

dG. 

B„ = 


dU 


(2-14) 


(2-15) 


Letting 


At 


r dU' 


\dt J 


= AU n 


(2-16) 


At 


r dH Y 


\dt J 


= AH n 


(2-17) 


( dUY 

At— = 8U 

\ dt J 


(2-18) 
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an implicit difference approximation to Equation (2-12) is 

. I + At —— + At 8£/” + ' = AU" +AH" (2-19) 

A z Ar 1 . 

where 

A = 4 

B = B,-B V 


In the above equation, D indicates a difference operator. The forward difference 
operator D+ and the backward difference operator D., are represented as follows: 



(2-20) 

DAA), J =U)„., -U),J 

(2-21) 


(2-22) 


(2-23) 


In order to more realistically approximate the physics of the governing 
equations, improve efficiency for the implicit scheme and reduce the numerical errors, a 
flux vector splitting technique is used in this study. Steger and Warming pointed out 
that 1411 

F = AU (2-24) 

and 


G = BU 


(2-25) 
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The Jacobian matrixes A and B con be diagonalized by similarity transformation S x 
and S y as follows 


\u 

>|° 

I o 

Lo 


0 0 

u + c 0 

0 u 

0 0 


0 1 

: k 

u-c\ 


(2-26) 


b = s;'w-' 


0 0 

V 0 

0 v + c 

0 0 


0 

0 


v-c 


with 


(2-27) 


1 0 0 

0 pc 0 

0 0 1 

0 - pc 0 


-1/c 2 

1 

0 

1 


(2-28) 



l 

0 

0 

-1/c 2 

s = 

0 

1 

0 

0 

y \ 

0 

0 

pc 

1 


_0 

0 

-pc 

1 


(2-29) 


1 0 0 0 
-uf p 1/p 0 0 

-v/p 0 1/p 0 

a(3 -wp -vP p 


(2-30) 


where c is the speed of sound, a=l/2(u 2 +v 2 ), P=y-1, and y is the ratio of the specific 
heats of the gases. In general some of the elements of the diagonalized matrix above are 
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positive and others are negative. Their signs determine the direction of information travel. 
So we can define the following matrices :* 26,271 


A = A^+A_ 

w+M 


=s;V" 


o 


o 


«+c+|w+c| 

0 r 1 0 


Saw 


1...-1 


0 
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2 
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0 

u + c — |w + c) 
2 

0 


2 
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2 
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0 

0 
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2 - 

0 

0 

0 

u-c— \u— c| 


IvuS. 


B = B^+B_ 

v+M 
2 

0 


=s;V 


s;V 


o 

o 

y-M 

2 

0 

0 
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0 

v+M 

2 

0 

0 

0 

v-M 

2 

0 

0 


0 

v+c+|v+c| 

- 

0 


0 

0 

0 


wS y + 


0 

V + C~M + c| 
2 

0 


v — c + |v — c| 
2 - 

0 

0 

0 

v— c— M~c| 


KmS„ 


(2-31) 


(2-32) 
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Using the direction of information travel, the flux crossing the surface separating 
volumes i,j and i+1 j is 



(2-33) 

= 5 + U, , + B.U^.j 

(2-34) 


Considering the flux vector splitting technique and performing the operations 
indicated in Equation (2-19 ), we obtain the following equation in terms of block matrices 
Cj, C2, C3, C4 and C5I 


W*. + C,SU,J + + CfiU MJ + CfiU MJ = A U’j + A HI (2-35) 


where 


C, =I + Ai\ 




( B_ B ' MN . MN . 
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VAr Ar Ar Ar 
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Equation (2-35) can be solved by line Gauss-Seidel iteration with alternating 
sweeps in the backward and forward z-directions, p9) as shown in Figure 2-1 . 


For k=l,3, ... 
Backward Sweep: 


C 2 iU‘j„ + C.St/.'j + CjiUfj., + C,Wl + C£U‘_l = At/”,- + AH-J (2-36) 


and 


For k=2, 4, ... 
Forward Sweep: 


C 2 st/,*j;, + Cfiutf + C,su‘j!, + CjbU&j + Cfiu^j -At/* + a//; 




(2-37) 


The above implicit block matrix equation can be put the following form: 
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where 


T 1 0 0 O' 

! 0 1 0 0 

F =1 

1 0 0 -1 0 ' 

Lo 0 0 1. 

The last equation in Equation (2-38) is used to set the boundary conditions 
which will be discussed in the later section. 

With the above numerical schemes the Upwind, Hybrid Explicit-Implicit 
equations have the following form: 

The predictor step 


A ur = -A/] 


M | P.(G,-G.) 


Az 


Ar 


-H 


<,j 


D..A D .A + A D + .B, 
I+At(-^ r ^ + ~^ Z ± ) + At(- - 


Az 
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Ar 
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ST/"; 1 = A U”j + AH?J 


J >7 


w+1 




(2-39) 


The corrector step 
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At/,";' = -a/1 


D.Jv £>.(<?, -G.) 


>+l 


Az 


Ar 


-// 


<.7 


D> £) .A D..B D .B. At T 
/ + A/(^ r - + -^) + A/(-^ + ^)--(8G v )J SC/”* = At/"J + AH”j 


A z Az 


u u =2 (U ‘-' +U % +W U> 


(2-40) 


where 5Ujj n+1 is calculated from Equation (2-38). 


2.3 Turbulence Model 


Turbulence modeling is the most important factor influencing the convergence of 
the Navier-Stokes equation solver, which is classified according to the number of 
supplementary partial differential equations. This number ranges from zero (algebraic 
model) to two (k-s model). The turbulent shear stresses in the mean-momentum 
equations are replaced by product of an effective viscosity and a mean rate of strain. The 
zero-equation model uses algebraic formulae to find the turbulent viscosity, which 
involve only properties of the mean velocity profile as unknowns. This is implied that the 
mean motion is unaffected by turbulence intensity and the length scale, which can be 
specified by an algebraic equation. One- and two-equation models obtain the velocity 
scale from a solution of the modeled form of turbulent kinetic energy equation and the 
specified length scale equation. In this study, the algebraic model has been used to yield 
faster convergence at reasonable accuracy. 

The algebraic turbulence model used in this analysis is a two-layer algebraic eddy 
viscosity model proposed by Baldwin and Lomax. 151 The effects of turbulence are 
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simulated in terms of the eddy viscosity coefficient, p.„ which is calculated for an inner 
and an outer region 

** ~\V-ouur Z > Z b 

where z is the normal distance from wall, and z ^ is the smallest value of z at which 
values from the inner and outer formulas are equal. 

For the inner region, the Prandtl-Van Driest formulation for turbulent viscosity 

is used: 

(P, );„„„ = pA: 2 /[l-e ( ^ ) ] 2 |©| (2-41) 


where co is the vorticity which is given by 


|co| = 


If du dvY 


and 



A + =26 


where x w , p w and are the local shear stress, density, and laminar viscosity 
evaluated at the wall. 

For the outer region, the Clauser formulation for turbulent viscosity is used: 


( Mv Pouter kC cp P ^ wake ^kleb O') 


( 2 - 42 ) 
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where k is the Clauser constant, C cp is an additional constant, and 


F wake = min 


f y max Fk 


MAX 


l ^WK y MAX H DIF I 


MAX - 


(2-43) 


The quantities y MAX and Fmax are determined from the function 

F(y) = >’|co|[l - e ( ^* ) ] (2-44) 


In wakes, the exponential term of Equation (2-44) is set equal to zero. The 
quantity Fmax is the maximum value of F(y) that occurs in a profile and y MAX is the 
value of y at which it occurs. The F k)eb (y) is the KlebanofF intermittency factor given by 


r M (y) = 



(2-45) 


The quantity u DIF is the difference between maximum and minimum total 
velocity in the profile 

u DIF = (Ju 2 + v 2 ) - Cju 2 +v 2 ) (2-46) 

\ ^max V. x m i n 

The constants used for this model have been determined by requiring agreement with 
the Cebeci formulation for a constant pressure boundary layer. The values determined are 
C cp =1.6, C wk =0.25, C kleb =0.3, k=0.4 and K=0.0168. 

In effect the voticity co is used to determine the length scale, so that the necessity for 
finding the outer edge of boundary layer is removed. This model has the advantage of 
avoiding the necessity for finding the edge of the boundary layer and exhibits good 


accuracy. 
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backward sweep direction forward sweep direction 


Figure 2-1 Sweep direction for line Gauss-Seider iteration (From Reference [27]). 



2.4 Diffusion Model for Thermal Radiation 


Under high temperature conditions of gaseous core reactor the radiative heat transfer 
rate in these opaque gases is comparable with the convective heat transfer rate. Therefore, 
the heat transfer analysis of a fissioning gas must include both convective and radiative 
heat transfer. The flow and radiation transport equation must be solved simultaneously to 
determine the temperature distribution and heat transfer rate. For such numerical model, 
radiative and convective heat fluxes are combined in the energy equation. Although there 
are several ways to obtain the radiative heat flux, for example, the Monte Carlo method, 
the zoning method, and the p-n approximation, the Rosseland diffusion approximation is 
widely used because it is simple to formulate and calculate. In order to use the Rosseland 
diffusion approximation the following assumptions are required: 

(1) The flow is an opaque and gray dense medium with only absorption and 
emission. 

(2) The radiation arriving at any location comes only from the immediate 
surroundings because any other radiation is absorbed before arriving at that location. 

(3) The particles are locally in thermal equilibrium and near velocity equilibrium. 

(4) The axial radiative dissipation is neglected. 

Under above assumption the radiative heat flux is proportional to the temperature 
gradient and can be written as:^ 47] 


?r 


= -K r S7T 


(2-47) 


[ 40 ] 


where K r is the radiative conductivity and is defined by 1 
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where a R is the Rosseland mean opacity and a is the Stefan-Boltzmann constant (5.67* 10' 
8 W/m 2 .K 4 ). To obtain the Rosseland mean opacity a R , the spectral absorption coefficient 
is averaged over the entire frequency range as the following:* 1 11 


) 


8B 

v 

6T 


dv 


<>r = 


1- 


dB„ 


a., dT 


dv 


(2-49) 


where v is the photon wave number, B v is the Planck function, a^, is the spectral 
absorption coefficient, and n is the real index of refraction. Direct calculation of the 
Rosseland mean opacity for UF 4 is difficult due to the lack of experimental data for the 
spectral absorption coefficient of UF 4 . However for this study, the Rosseland mean 
opacity can be estimated by* 12 ^ 


°r = n ° p h ( 2 ’ 5 °) 

where N is the molecular number density of the gas and a ph is the photon collision cross 
section per molecule. For UF 4 gas a ph is estimated to be equal to 2.76x10' m . Using 
the assumption of the perfect gas law the molecular number density can be estimated as 

W = 4 (2-50 


where k is the Boltzmann constant (1.3806*10' 23 J/K), p is the gas pressure, and T is the 
gas temperature. Therefore, the radiative conductivity can be written as 


16 ctk , 
K r = ^Z — T* 

^ P hP 


(2-52) 


Using the estimated opacities the Rosseland approximation leads to a considerable 
simplification in the expression for radiative heat flux. The radiative heat flux is added to 
the Equations (2-3) and (2-4) to establish the balance of the energy equations. Since 
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radiation leaving from the surface is not taken into account, the coefficient of radiation 
conductivity in Equation (2-52) has to be changed from 16/3 to 8/3 


2.5 Wall Heat Transfer Model 


The system is limited to smooth, straight tubes within which the fluid flow is 
turbulent. In the system, two boundary conditions, constant wall temperature and constant 
wall heat flux on the outside surface of the tube, are of particular interest and considered 
here. There are two ways to describe the convective heat transfer into the wall boundary. 
One way is to use Newton’s cooling law. That is 

r = h(T,-T t ) (2-53) 


where h is the convective heat transfer coefficient or so-called unit thermal conductance, 
which is not a material property as thermal conductivity. It is a complex function of the 
composition of the fluid, the geometry of solid wall and the hydrodynamics of the fluid 
motion, in particular, the temperature distribution near the solid wall. So it is impossible 
to estimate the heat transfer coefficient without solving the Navier-Stokes equations. 

The other way is to apply Fourier’s law. That is 


q" 


dT 

v \dr ) . 


gas 


(2-54) 


where k w is the thermal conductivity which is evaluated at the wall temperature. The 
temperature gradient in equation (2-54) is calculated using the wall temperature and 
adjacent cell temperature. It may be numerically written as 


q n = -k y 


T - T 

£2 1 w 

r 2 -r. 


(2-55) 
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The condition of the above equation is that the temperature distribution near the wall 
has to be linear. This means that equation (2-55) can only be applied in the viscous 
sublayer. 

The turbulent boundary layer is composed of three different sublayers, such as 
viscous sublayer, buffer layer and turbulent core, which is presented in Figure 2-2. In the 
viscous sublayer, heat transfer is dominated by diffusion. The temperature distribution is 
linear, as shown in Figure 2-3. In order to use equation (2-55), it is important to 
determine the thickness of the viscous sublayer which may be estimated by the following 
relation: 135 ^ 


72.94z 



(2-56) 


The magnitude of 8 b is about 2-5 pm. This means that at least 3 cells have to locate in 
the 2-5 pm region near the solid-fluid interface. A simple sketch of finite volume mesh to 
express the boundary index system is shown in Figure 2-4. 


2.6 Boundary Conditions 


Following the above wall heat transfer model, three kinds of thermal boundary condition 
are applied in this study, which are adiabatic boundary condition, a constant wall 
temperature condition and constant heat flux condition. 

The adiabatic boundary condition 

The general form and numerical expression are 
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q” = 0 

(2-57) 

t, = t 2 = t w 

The physical meaning of this boundary condition is that no heat transfer passes through 
the boundary. 

The constant wall temperature condition 
The general form and numerical expression are 

T v = f{z) (2-58) 

T x = 2T W - T 2 (2-59) 

The physical meaning of this boundary is that temperature is on the isothermal 
condition. This boundary condition can be specified for the solid-fluid interface, which is 
indicated in Figure 2-5. 

The constant heat flux condition 

The general form and numerical expression are 


q ” = const 


(2-60) 


T x 


- T 
1 2 


+ 


li 

k 


(r 2 ~r x ) 


(2-61) 


Because both T[ and T 2 are unknown before each iteration, it is difficult to specify the 
boundary condition which is indicated in Figure 2-6. An enthalpy-rebalancing scheme is 
developed to apply this boundary condition. It is performed in the next section. 
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2.7 Enthalpy-Rebalancing Scheme 

For the problem with constant heat flux condition, which is shown in Figure 2-4, it is 
difficult to obtain convergent solutions. The reason is that integrating solution is 
impossibly found just from its first derivative. In order to solve problems under an 
arbitrary heat flux boundary condition, a novel method of enthalpy-rebalancing at the 
transverse flow surface has been developed. This scheme is based on the fact that under 
steady conditions the gas enthalpy rise at each transverse flow surface is equal to the heat 
removal from the wall. Mathematically, it can be written as follows: 

(AQ)' = 2nRAzq" (2-62) 

Equation (2-62) is explicitly solved to obtain the balancing bulk temperature (T b )j , 
that is 



(2-63) 


where T b is the bulk temperature, and Cp, is the local specific heat. In this section, (T b )j 
is known. (T b ) i+1 can be found from Equation (2-63) which is always true in the transient 
procedure based on the enthalpy-balancing principle. The bulk temperature ((T b ) i+1 ) num 
is calculated from temperature fields under every iteration. Since these two bulk 
temperatures should be equal, their relation may be used as a convergence criterion as 
follows: 


6 = 


(5k 

((U + i)™ 


(2-64) 



32 


Before achieving the steady-state flow conditions, the mass flux pu, is always smaller 
than the final steady value. Therefore, during the transitory pre-steady state conditions the 
value of 5 is greater than one. Whenever the mass flux approaches its steady-state value, 
the convergence parameter 5 approaches its asymptotic value which is one. At this 
moment, the wall boundary condition (2-61) is used to calculate the wadi temperature. As 
soon as 6 is equal to one, a thermal steady-state condition is achieved. At this point, the 
wall temperature is fixed to control the heat transfer through the solid wall until a global 
thermal convergence is achieved. This scheme has been proven successful when 
problems with heat flux boundary condition are to be solved. 


The temperature distribution in the fuel elements is the essential to the prediction of 
the lifetime behavior of these components. The temperature distribution in the fuel 
dependents on the volumetric heat generation rate, fuel material properties, dimension of 
the fuel element and geometric configuration. Of these the first 3 parameters are related 
the neutronic design considerations, so geometric configuration is the only optional factor 
via a thermal analysis. There are four types of geometric configuration, plate, pin, tubular 
and square lattice honeycomb, which are commonly used in reactor fuel element design. 
The configurations can be found in Figure 2-7. In order to evaluate the geometric 
configuration and predict the temperature distribution, a dimensionless number is defined 
as the following: 


6(a) = 


kAT 

q m D h ( a ) 2 


(2-65) 
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where q'" is the volumetric heat generation rate, D h is the hydraulic diameter, k is 
conductivity of the fuel materiel and AT is the difference between maximum temperature 
and wall temperature and a is void parameter which is defined as 

^Fiov ( 2 - 66 ) 

A 

" Tola f 

The physical concept behind this definition is the efficiency of the conductive heat 

transfer. The focus of the Delta number is on the geometric behavior of the fuel element 
which is strongly dependent on the temperature distribution. The thermal conduction 

equation describes the temperature distribution in the fuel element. Under the assumption 
of negligible thermal expansion, the general equation of heat conduction becomes : [3 11 

V 2 T(x,y)+^- = 0 (2-67) 

where T is temperature distribution function (K) 
q'" is the volumetric heat generation (w/cc) 
k is thermal conductivity (J/cm/K) 

2.8.1 Plate and Pin Configurations 

Because of the symmetric geometry, at steady state the general conduction equation 
for geometries of plate and pin reduce to a ordinary differential equation which can be 
written as: 

For plate. 
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For pin 


(2-68) 


1 d ( dT(r)) 

r 

r dr v dr J 




(2-69) 


Equations 2-68 and 2-69 can be integrated, respectively: 

rt*)=r„+£(/ : - 4 * ! ) 

^ )=r - + S ( ° ! - 4r2) - 

Therefore, 

. d-q ) 2 

na,e 32a 2 

. 0-a) 2 

Pm 16a 2 

2.8.2 Tubular and Square Lattice Honeycomb Configurations 


(2-70) 


(2-71) 


(2-72) 


(2-73) 


For complex geometry, such as tubular and square lattice honeycomb, partial 
differential equation (2-55) can not be solved analytically. The computer code AN SYS is 
used to get the numerical results. 
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The general purpose Finite Element code ANSYS developed by Del Salvo et al. uses 
a set of elements whose stiffness matrices have already been set up. These elements are 
referred to by a pair of numbers, each identifying the type of analysis and the degrees 
of freedom per node. ANSYS uses the following relationships for heat transfer analysis: 

V 2 r(xj>)+y = 0 (2-74) 

These equations are solved by the Finite Element (FE) technique. To start with, a 
temperature distribution is obtained based on an input reference temperature and then 
properties such as thermal conductivity, are evaluated at nodal temperatures. This 
procedure is continued until a convergence in nodal temperature is achieved. Physical 
properties of the structural material at each node are then evaluated based on this 
calculated nodal temperature. These nodal temperatures are stored in a file for later use. 

The input parameters and sample input file are listed below: 

Table 2-1 Input parameters for ANSYS analysis. 


Parameter 

Symbol 

Value and Unit 

Volumetric heat-generation rate 

q’" 

9 (W/mm 2 ) 

Surrounding temperature 

T s 

2500 (K) 

Heat-transfer coefficient 

h 

0.0655 (w/mm 2 .K) 

Thermal conductivity 

k 

0.036 (W/mm.K) 

Specific heat, constant pressure 

c P 

4.44 (J/kg.K) 



















The sample data for ANSYS analysis: 


/title, tubular 

kan, -1 

treff,2500 

et,l,55 

iter,- 1 

ktemp,-l 

k, 1,0.0 

k,2, 0.278 

k,3,0.62 

k,4, 0.791 

k,5 ,0.791, 0.53 

k, 6, 0.371, 0.6425 

k,7,0. 1855,0.321 3 

k,8,0. 11,0.19 

L,l,2 

L,2,3 

L,3,4 

L,4,5 

larc,5,6,4,-0.84 

elsiz, ,50 

L,6,7 

L,7,8 

L,8,l 

a,l,4,5,6 

amesh,all 

wsort,x 

wsort,y 

lssel, line, 1,4, 1,1 
hflow, all, heat, 0.0 
nail 
lsall 

lssel, line, 5,„1 
eall 

cvsf,all,z„0.0655,2500 

nail 

lsall 

lssel,line,6,8,l,l 

hflow,all,heat,0.0 
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nail 

eall 

kxx,l, 0.036 

kyy, 1,0.036 

c, 1,4 .44 

q,all,9 

Iwrite 

afwrite 

FINI 

The Figures 2-8 to 2-12 show the typical computation grids and temperature distribution 
for tubular and square lattice honeycomb, respectively. The numerical results are 
presented in the following tables: 

Table 2-2 Variation of Delta number in tubular fuel element. 


a 

D h 

AT(K) 

5 

0.1 

1.021 

162 

0.7826 

0.2 

1.286 

111 

0.2681 

0.3 

1.576 

72 

0.1161 

0.4 

1.822 

49 

0.0592 

0.5 

2.034 

33 

0.0319 

0.6 

2.228 

24 

0.0193 

0.7 

2.407 

17 

0.0117 

0.8 

2.573 

14 

0.0085 
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Table 2- 3 Variation of Delta number in SLHC fuel element. 


a 

D h 

AT(K) 

6 

0.1 

0.947 

298 

1.3297 

0.2 

1.341 

188 

0.4178 

0.3 

1.643 

128 

0.1896 

0.4 

1.897 

88 

0.0977 

0.5 

2.121 

59 

0.0525 

0.6 

2.324 

38 

0.0282 

0.7 

2.512 

22 

0.0141 

0.8 

2.683 

9 

0.0052 


After polynomial fitting, the relations between Delta number and void parameter are 
given as follows: 


Table 2-4 Polynomial formula of the Delta number. 


Configuration 

Delta Number 

Plate 

(1-a ) 2 


32a 2 

Pin 

(1-a) 2 


16a 2 

Tubular 

2.1 - 1 8.9a + 73.6a 2 - 142.1a 3 + 135.7a 4 - 50.6a 5 

SLHC 

3.8 - 36.5a + 147.2a 2 - 294.4a 3 + 286.7a 4 - 109.4a 5 


The above relations are presented in Figure 2-13. 
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2 - 8 -3 Applicatio ns of Delta Nnmhpr 

The delta number is similar to other dimensionless numbers in heat transfer, such as 
Mach number and Reynolds number. It expresses the relationship between several 
thermal parameters which are heat-generation rate, thermal conductivity, surrounding 
temperature, maximum temperature and the geometric parameter. There are at lest two 
important applications related to design of the nuclear fuel element: (1) to evaluate the 
coolant channel configuration and (2) to predict the maximum temperature. 

Delta number is function of the void parameter for several different coolant channel 
configurations. From conduction of view, the optimum configuration can be determined 
based on Figure 2-8. The smaller Delta number means the lower maximum temperature 
in the fuel element. So the Delta number can be used to control temperature gradients, 
which is important for the lifetime of reactor components. 

Once the fuel surface temperature and volumetric heat generation rate have been 
determined, the centerline temperature, maximum temperature of fuel element, can be 
calculated from equation (2-65). This method to calculate centerline temperature, instead 
of the time-consuming CFD and other iterative procedures, will remarkably simplify the 
centerline temperature calculation. This method has been used by Reference [15] and 
proved to be accurate and efficient. 




Turbulent Core 


Butler Layer 


Viscous Sublayer 


Figure 2-2 Boundary layer of turbulent How (From Reference (42]). 



Fiimre 2-3 Turbulent boundary layer and laminar sublayer. 
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entrance 


q"(z) 



exit 


Figure 2-6 The heat flux boundary condition. 
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(1) Plate 


(2) Pin 


(3) Tubular 


(4) Square Lattice 
Honeycomb 


Figure 2-7 Configuration of reactor fuel elements. 
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Figure 2-8 Temperature distribution in square lattice honeycomb fuel element (a is 
ranged from 0.1 to 0.8). 
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Figure 2-10 Temperature distribution in square lattice honeycomb fuel element (a=0.2). 
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Figure 2-1 1 Temperature distribution in tubular fuel element (a=0.2). 



48 



Figure 2-12 Temperature distribution in tubular fuel element (a=0.2). 



Delta Number 


49 






CHAPTER 3 

THERMOPHYSICAL PROPERTIES 


Thermophysical properties are very important parameters in CFD calculation. In 
order to get accurate results, the real thermophysical properties of hydrogen, helium and 
UF 4 gases are used in this study. 


3.1 Dissociated Hydrogen 

The data of the hydrogen property package comes from three different sources which 
provide information in different temperature regions. From temperatures ranging from 
13.8° to 3000 °K, data is based on the National Bureau of Standards. Above 3000°K and 
up to 10,000°, the data is obtained from NASA computer programs. Above 10,000°K, the 
data is generated by extrapolation of the data base. For the hydrogen properties, the data 
provided by the hydrogen package is the best available at this time. 

In order to get the hydrogen properties for given condition, a numerical code based 
on cubic spline interpolation was developed. This code interpolates with respect to 
temperatures or enthalpies and pressures to find all other thermodynamic properties, such 
as density, viscosity, thermal conductivity, entropy, specific heat, speed of sound and 
ratio of specific heats. 

The concept of the cubic spline interpolation is to construct a cubic function S k (x) on 
each interval [x k , x k+1 ] so that the resulting curve y = S(x) and its first and second 
derivatives are all continuous on the large interval [x^, x N ]. The function S(x) is called a 
cubic spline which has to satisfy the following equations: 1361 
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I. S(x) = A k +B k (x-x k ) + C k (x-x k ) 2 + D k (x-x k ) 3 

(3-1) 

II. S(x) = y k 

(3-2) 

III .S k (x M ) = S M (x M ) 

(3-3) 

iv. su** + i)=s; + ,(x* +1 ) 

(3-4) 

v. 5' k (x jl+1 ) = £*+i(**+i) 

(3-5) 


The above properties of cubic spline function are shown in Figure 3-1 . 


The general form of cubic spline equation S k (x) is the following expression: 


S k (x) = 


S"(x k )(x k+l - x) 3 S'(x M )(x t - x) 3 


6(** +1 ~x k ) 


6(** + , - x k ) 


+ (- 


y k 


S”(x k )(x k+l — x k ) 


** + i -*k 


)(x k+] -x)~ 


T* + i S”(x k+l )(**♦! -x k ) 


X M -x k 


)(x k - x) 


(3.6) 


Now use the properties IV and V to obtain an important relation with respect to x k , 
x k+ i,y k and y k+1: 


S%Xji - 1 )(*& x k _] ) S\Xji )(XjUi x*_i) & (x_k+i )(**+i x - k ) 

6 + 3 6 
(y™ zzA ( ^ -y*-i) 

** +1 -x k x k -x k . t 


(3.7) 


This equation is a system of n-1 linear algebraic equations to be solved simultaneously. 
Its coefficient matrix is tridiagonal matrix. The linear system is diagonally dominant and 
has a unique solution. After the linear system is solved, the unknowns, S”(x k .,) , S”(x k ) 
and S”(x k+1 ), are determined. The interpolation value can be calculated from equation 
(3.7). 

This code was used to draw 8 sets of curves that look very smooth when viewed by 
the eye. The curves are shown in Figures 3-2 to 3-9. 
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3.2 Helium Gas 


The thermophysical properties, heat capacity C p , viscosity p and thermal 
conductivity IQ, based on the references [34] and [45] are used throughout the 
calculation. The formula of polynomial fitting are generated as follows: 


C(kJ / kg. K) = -8.3 799x1 O' 8 T 2 +351 145x1 O' 4 T + 0.87076 


\i(N.s / m 2 ')= 


3.5761x10 T - 8.4466x10 7 1 +0.909 139T + 1.59283 

1.0 7 


K C (W / m.K) = 


0.257455T + 95.4545 
10 3 


(3-8) 

(3-9) 

(3-10) 


3.3 Uranium Tetrafluoride Gas 

Thermophysical properties of UF 4 gas used throughout the calculation are given 
by Equations (3-ll)-(3-13). [3] In these equations, C p is heat capacity, p is dynamic 
viscosity and K c is thermal conductivity. 

Analysis of existing data and rigorous theoretical calculations are used to developed 
the following heat capacity equations for the gaseous UF 4 in 1000° to 10,000°K range. 

C p (J / mol.K)= 1215 + 2.24 x 10" 3 T-— [T<3500°K] (3-11) 


C p (J / mol. K) = 124.12 - 1.28 x 1 0' 3 T 


[ T > 3500°K ] 


(3-12) 
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Using the semi-empirical relations for related transport parameters, the viscosity of 
UF 4 in 1000° to 10,000°K range is developed [30] as follows: 




3.357 x 10 -6 Vr 

a + bT 


(3-13) 


where 


a = 0.8 b = -7.1x10' 


[ T < 3500°K ] 


a = 0.67 b = -2.04 x 1 0' s [ T > 3500°K ] 


Thermal conductivity of UF 4 gas is estimated as: 
K C (W lm.K)= 3.2 p (C p + 1 0.393) 


(3-14) 


where heat capacity C p is in J/mol.K, and the viscosity p is in Pa.s. 
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Figure 3-7 Specific heat of dissociated hydrogen for different pressure. 











CHAPTER 4 

ASSESSMENT OF NUMERICAL MODELING 


In order to check the assessment and validation, several problems have been carried 
out to compare with published results. Excellent agreement with theoretical and 
experimental results has been obtained. The details of the tests are indicated in the 
following sections. 



Fluid flow studies over a sphere which is bounded in a cylindrical tube is of 
fundamental importance and has received attention for bounded flow. This is typical 
model for the low-drag measurement and determination of fluid viscosity. One of the 
parameters of paramount importance for the fluid flow over a bluff body is the separation 
angle. It is required for the drag force calculation and very difficult to ascertain by 
experimentation. In this study a numerical analysis is performed to observe the 
dependence of the separation angle on the diameter ratio between the sphere and the 
cylinder (y). This study is based on the assumption that the Reynolds number is fairly 
high but low enough to have a laminar boundary layer. Figure 4-1 to 4-4 illustrate the 
flow patterns and separation angles for the cases y=0.3, y=0.4, y=0.6 and y=0.8, 
respectively. The results show that the separation angle reaches a minimum around y=0.5. 
This result provides very good agreement with following semi-empirical equation which 
is based on the boundary layer analysis: 1351 
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<£(0.9277) + <£ 3 (2.095/ 2 -0.364l)+<£ 5 (3.2466/ 4 -15495/ 2 +0.039)+ 

+<£ 7 (4.396/ 6 - 1 3.4907y 4 + 0.4687y 2 - 0.0015)= 0 ' 5 1 * 

where <£ is function of separation angle. The following table provides the numerical 
results to compare the boundary layer approximation: 


Table 4-1 Comparison of separation angle between numerical calculation and boundary 
layer approximation 


DIAMETER RATIO 

SEPARATION ANGLE (DEGREE) 

y 

Numerical Calculation 

Boundary Layer Approximation 

0.3 

110.0 

108.46 

0.35 

107.8 

102.59 

0.4 

105.8 

92.61 

0.45 

103.9 

90.33 

0.5 

104.0 

88.35 

0.55 

105.7 

88.98 

0.6 

107.6 

90.53 

0.65 

109.1 

96.03 

0.7 

110.5 

106.16 


The comparison is shown in Figure 4-5. 

An experiment was performed to measure drag force on a ball by using a Cahn 
instrument manufactured C2000 model microbalance. Teflon balls of different diameter 
were centrally hung in a long cylindrical tube of diameter 22.54 mm. The spheres were 
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suspended by means of a fine nylon string (d=0.07 mm), which can be found in reference 
[24]. A numerical investigation was carried out to compare with the experimental results. 
The drag force can be found from the following formula: 



where p.; is viscosity. 

Figure 4-6 provides the comparison between the numerical and experimental results. 
It is shown that by increasing the number of grid points on the sphere, the drag force 
results finally approach the experimental values. 


Among the long list of correlations for the Nusselt number, axial distance corrections 
and property corrections, the following correlations which have the more relevance to 
heat transfer in a nuclear reactors are chosen to justify the validity of computational 
results.^ 15,38 ' 


(1) Empirical Correlations of Nusselt Number 
1 . Colburn correlation (1933) 
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6. Notter-Sleicher correlation (1975) 
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7. Churchill correlation (1977) 
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(2) Axial Distance Corrections for Nusselt Number 


1. Perkins and Worsoe-Schmidt correction (1965) 
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2. Pierce correction (1963) 


f : = 1 + 0.3^ 


3. Bussard correction (1965) 
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4. Reynolds-Swearington correction (1965) 
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(3) Property Corrections for Nusselt Number 


1 . Notter-Sleicher property correction 
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3. Perkins and Worsoe-Schmidt property correction 
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(4-16) 


Figure 4-7 compares the numerical results with four Nusselt number correlations. 
These correlations include the Colburn equation, the Gnielinski correlation, the Karman- 
Boelter-Martinelli equation and the Notter-Sleicher correlation. Figure 4-7 shows that 
numerically calculated Nusselt number with some of empirical correlation results is 
almost undistinguished. This comparison not only evaluated the empirical correlations of 
Nusselt number, but also assessed the validation of numerical modeling. 






Figure 4-4 Flow pattern and separation angle for y=0.8 (0=1 14.3°) 
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Figure 4-5 Comparison of separation angle between numerical calculation and boundary 
layer approximation. 
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Figure 4-7 Comparison of the numerically calculated Nusselt number with the 
empirical correlations. 




CHAPTER 5 

RESULTS AND DISCUSSION 


This chapter is composed of three parts. The first part presents the numerical results 
for the fluid flow and heat transfer in the core of a gaseous core nuclear reactor. UF 4 gas 
is used to simulate the flow of a Ultrahigh Temperature Vapor Core Reactor with 
Magnetohydrodynamic generator (UTVR-MHD) system. The second part describes the 
numerical results for fluid flow and heat transfer in the core of a gas cooled reactor. 
Helium gas is utilized in this part. The last part indicates the results for the fluid flow and 
heat transfer in XNR2000 nuclear rocket core. Hydrogen gas is considered in this study. 


A gas core reactor uses a gaseous fissile material as fuel for power generation. This 
allows for power generation at temperatures much higher than the melting point of solid 
fuel nuclear reactors. Power generation and power conversion at very high temperatures 
can potentially reduce the system mass and improve the specific impulse of a nuclear 
thermal rocket. One of the most challenging issues related to design and operation of gas 
core reactors is the containment of the fissioning plasma. The wall cooling is the most 
important issue in design of gaseous core nuclear reactors. The heat transfer process 
involved in an ultrahigh temperature gas core reactor systems is characterized by the 
convective flow of a radiating gas. Uranium and uranium compound gases even at pressures 
in the range of 10 to 40 atm are optically thick. At temperatures close to 3500°K which is 
the typical exit temperature of the reactor core in a recent design, the radiative heat 
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transfer rate in these opaque gases is higher and comparable with the convective heat 
transfer rate. 1121 Therefore, the heat transfer analysis of a fissioning gas must include both 
convective and radiative transfer. The flow and radiation transport equations must be solved 
simultaneously in order to determine the temperature distribution and heat transfer rate. 139,421 

The model presented in this paper considers combined convection and radiation of 
compressible, turbulent and developing flow of a radiating gray gas under very high 
temperature conditions. Although the combined convective and radiative heat transfer 
problem has been the target of many investigations, a comprehensive numerical solution of 
this problem has not yet been published. 121 ’ 321 

The case selected for this analysis consists of a cylindrical tube lm long and 0.05m 
in diameter, i.e. z/D=20. The constant wall temperature is set at 1600°K, and total inlet 
temperature is specified as 2000°K, which is a typical design temperature at the core inlet 
of the ultrahigh temperature vapor core reactor system. A stagnation inlet pressure of 2 
MPa, the back pressure of 1 .8 MPa and three different values of internal heat generation 
rate of 100 MW/m 3 , 500 MW/m 3 and 1000 MW/m 3 are used. The Reynolds number is 
ranged from 10 4 to 10 6 . 

The flow domain is divided into 54 radial and 54 axial control volumes spaced non- 
uniformly. A fine grid, defined by algebraic method, was used near the wall to ensure 
more than two grid points in the laminar sublayer, as shown in Figure 5-1 . 

5.1.1 The Temperature Distribution 

The temperature distributions, which are obtained at a constant wall temperature of 
1 600°K and inlet stagnation temperature of 2000°K, are used for verification of effect of 
internal heat generation rate and are presented in Figure 5-2. The relations between 
temperature distribution and internal heat generation rate will be discussed in the next 
section. The numerical results are compared with the results of Hoogenboom et al. 1 ’ 81 
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Figure 5-3 shows that the maximum temperature and the trend are almost 
undistinguished. The Nusselt number for pure convection is compared with the results of 
Petukhov-Kirillov-Grielinski equation, which is shown in Figure 5-4. The differences 
between those two results are less than 3% at the tube exit (z=lm). The excellent 
agreement between the present model and the previous literature indicates the validity of 
the numerical modeling for gaseous core reactors. 

All of the above results are obtained from the solution of the flows at the steady- 
state. In order to verify the steady state solution, the convergence history for the velocity 
and temperature residual is examined, as shown in Figure 5-5. In this case, the steady 
state is reached when the root-mean-squared residuals of velocity and temperature drop 5 
and 4 orders of magnitude in about 900 iterations or 300 minutes of CPU time using a 
486-66 computer. The velocity and temperature residuals are defined by 


5 „ = 



(/maxX/max) 



(5-1) 


5 r = 



(imax)(/max) 


(5-2) 


5.1.2 Effect of Heat Generation Rate 


Gas core nuclear reactors are characterized by high internal heat generation rate in 
the fissile gas. The numerical model is used to predict the temperature distribution and heat 
transfer rate for the gases with internal heat generation. Uniform heat generation rates 
ranging from 100MW/m 3 to 1000MW/m 3 are included in the energy equation as source 


terms. 
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In Figure 5-6, the temperature distribution is presented for the three different values of 
heat generation rate. Because of the radiative heat transfer, the curves become more flat and 
grow faster than those of pure convection. The convective heat flux at the wall increases as 
the power generation rate increases. The Figure 5-6 also shows that the temperature at a 
heat generation rate of 1000 MW/m 3 increases more rapidly than others because the 
radiation at high temperature is more dominant. 

The maximum temperature, corresponding to the heat generation rates, are 2150°K , 
2750°K and 3550°K, respectively. This is indicates that a heat generation rate higher than 
1000MW/m 3 is necessary to maintain the gas temperature at about 3500°K, which is 
typical temperature required to achieve high efficiency in the gas core reactors. 1321 The 
bulk temperature gradient at the wall is steeper than it would be in pure convection, which 
is shown in Figure 5-2. These results indicate that the maximum temperature in the reactor 
is strongly dependent on the value of heat generation rate. 

5.1.3 Convective and Radiative Heat Fluxes 


The heat fluxes are presented for both convective and radiative heat transfer. The 
convective, radiative and total heat flux transferred to the wall are calculated using the 
following equations: 1371 
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where q t , q c and q r are the total, convective and radiative heat fluxes, respectively. The 
convective and radiative heat fluxes for the three different cases are presented in Figure 5- 
7. As the heat generation rate increases, the corresponding heat flux increases. It is found 
that at temperature close to 3500°K the radiative heat transfer rate is comparable with the 
convective heat transfer. It is seen that the radiative fluxes increase as higher wall 
temperatures. This behavior is due to the fact that the radiative wall heat flux is dependent 
on the local conditions at the wall. 


5.1.4 Nusselt number calculation 


Nusselt number is the parameter which characterizes the heat transfer rate at the wall 
boundary. The total, convective and radiative Nusselt number can be defined as: 


k(r, - r») 

q r D 


Nu = 


Nu = 


Nu , = Nu c + Nu r 


(5-6) 
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(5-8) 


Figure 5-8 illustrates the relation between local Nusselt number and Reynolds 
number. It is found that the Nusselt number increases as the heat generation rate 
decreases. This is mainly due to the increase of dynamic viscosity of UF 4 gas as 
temperature increases. Figure 5-8 also shows that the entrance length increases as the 


Nusselt number decreases. 
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Figure 5-1 Overall view of grid system. 
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Figure 5-3 Temperature distribution comparison between results obtained with the 
present methodology and by Hoogenboom et al. 
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Figure 5-4 Convective Nusselt number comparison between results obtained with the 
present methodology and by using Petukhov-Kinllov-Gnielinski correlation. 





esicliia 


xmperature-RBsidual 





• — Q”'= 1000 MW/m A 3 
f — Q'"= 500 MW/m A 3 
f — Q'"= 100 MW/m A 3 




Convective Heat Flux (MW/m 


86 



0 0.2 0.4 0.6 0.8 1 

x (m) 


Figure 5-7 Convective and radiative heat flux as a function of axial location: -x- Q= 100 
W/m 3 , -A- Q=500 MW/m 3 , -0-Q= 1000 MW/m 3 . 
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Figure 5-8 Local Reynolds number and Nusselt number variation: -x- Q= 100 W/m 3 , 
-A- Q=500 MW/m 3 , -Q-Q=l 000 MW/m 3 . 
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5.2 Calculations for a Gas Cooled Reactor 

Over the past few decades, high temperature gas cooled reactors (HTGRs) have been 
considered for a wide range of applications. Compactness, high efficiency and very high 
temperature capability of these reactors are of great importance to power generation in 
space. Other applications of HTGRs include: terrestrial electric power generation, nuclear 
thermal propulsion and direct use of high temperature gas for a variety of industrial 

MQ1 

processes such as steel-making. 1 The technology of HTGRs for commercial power 
generation, as well as the computational methods for analysis of thermal fluid 
performance of these reactors are well developed. The majority of existing HTGR 
thermal fluid analysis methods use empirical correlation to resolve heat and momentum 
transfer at the fuel-coolant boundary. However for HTGR concepts with operating 
parameters beyond those of commercial HTGRs, the issue of the accuracy and 
applicability of empirical correlations are not fully resolved. In particular, energy 
transport in very compact space power reactor concepts may require flow at very high 
velocities and high Reynolds numbers. This study presents a non-correlation based 
computational thermal-fluid model to analyze flow and heat transfer in HTGR cores. The 
computational model is also used to assess the performance of several mechanistic 
correlations for calculation of pressure drop in HTGR 

A dual path cermet fuel fast spectrum reactor is used as a computational model to 
analyze the high-temperature gas-cooled reactor system. The core consists of 631 fuel 
rods and 37 holes per fuel rod which has a flow-equivalent diameter of 3.2 mm and 
heated length of 0.544 m. The stagnation inlet pressure of 6.5 MPa, the back pressure of 5 
MPa and the core power of 25 MW, 50 MW, 75 MW and 100 MW are used, 
respectively. The Reynolds number of the flow is in the range of 10 4 to 10 5 . 

The flow domain was divided into 54 radial and 54 axial control volumes spaced 
non-uniformly. A fine grid, defined by an algebraic method, was used near the wall to 
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ensure more than two grid points in the laminar sublayer which is used to control the heat 
flux near the wall. Figure 5-9 shows the flow pattern which presents the typical velocity 
distribution of turbulent flow. 

5.2.1 Temperature Distribution 

The maximum temperature in the HTGR reactor is strongly dependent on the power 
density. To demonstrate the influence of power density, a calculation was done with core 
power of 25 MW, 65 MW and 100 MW. The temperature distributions under three 
different core powers axe shown in Figure 5-10. The maximum temperatures, 
corresponding with the core power, are 1500°K , 2200°K and 3000°K, respectively. This is 
indicates that a core power higher than 100MW is necessary to maintain the gas 
temperature at about 3000°K. Because of high heat flux the derivative of temperature at 
the wall is much higher than the derivative of temperature at the center area. The 
temperature at Q=100 MW grows much faster than temperature at Q=65 MW. The bulk 
temperature and pressure distributions for different mass fluxes are shown in Figures 5-1 1 
and 5-12, respectively. 

5.2.2 Nusselt Number Calculations 

Nusselt number correlations are of great importance to calculation of heat transfer. 
However, almost all of these correlations are developed under fully developed and 
constant wall heat flux conditions. In some cases the heat flux used for the generation of 
the experimental data base is rather low. The low wall heat flux indicates small 
temperature gradient in the flow boundary layer. Therefore the changes in flow 
properties due to temperature gradient may not be important. Under flow, temperature 
and heat flux conditions of ultrahigh temperature and compact HTGRs, similar to those 
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proposed for space power and propulsion applications, the flow is not fully developed, so 
the temperature gradient could be very large. So, the wall thermal boundary conditions 
may be different from what is used for the development of the experimental data base 
used for the derivation of the particular heat transfer empirical correlation. Thus, the 
detailed computational analysis developed in this work is used to evaluate the most 
common experimental correlations for wall heat transfer in HTGRs. 

Nusselt number variation along the tube length is shown in Figures 5-13 and 5-14. 
From the results of this study, the Nusselt number increases as the core power increases. 
The differences of the Nusselt number of core power of 25 MW, 50 MW, 75 MW and 
100 MW at x=0.554 m are about 14% and differences of the Nusselt number of mass flux 
of 45, 60, 75 and 90 kg/m A 2.s is about 29.5%. It is interesting to note that the core power 
and mass flux have similar effects to Nusselt number. 


5.2.3 Pressure Drop Calculations 


Similarly, simplified equations are used to calculate the pressure drop in heated 
channels of HTGRs. These equations may or may not be applicable to the flow and 
thermal conditions of ultrahigh-temperature gas-cooled reactors. An analysis is performed 
to evaluate the accuracy and validity of the pressure drop equations for HTGRs. 
Equations which are conveniently used to calculate the accelerational and frictional 
pressure drop in HTGR cores are: 
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Figure 5-9 Flow pattern for gas cooled reactors. 



Figure 5-10 Comparison of temperature distribution for different valued of core power: 
(1) Q=25 MW, (2) Q=65 MW, (3) Q=100 MW. 
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(5-11) 


Figure 5-15 shows the axial pressure distribution in the heated flow channel as 
calculated by the detailed Navier-Stokes solver and also by Equations (5-9), (5-10) and 
(5-11). In Figure 5-15, Incm implies equation (5-9) which is found in Reference [2]; 
Coml implies equation (5-10) and Com2 implies equation (5-11) which can be found in 
Reference [22]. From this analysis it is evident that the equation (5-9) provided the best 
agreement with the numerical results. This is mainly due to the fact that the thermal 
evolution of the flow is taken into account in the equation (5-9), which is much closer to 
the governing equation used for the CFD calculation. 


The XNR2000 is a fast-spectrum cermet-fueled nuclear reactor that heats hydrogen 
in a two-pass folded-flow configuration and delivers the hydrogen propellant to the 
nozzle chamber before expansion through a nozzle. 1141 The total engine flow rate is 12.1 
kg/s and cycle operates at a chamber pressure and temperature of 5.5 MPa and 2850°K. 
The engine delivers 111.2kN of thrust at a thrust to weight ratio of 5.1 and a specific 
impulse of 944 seconds including kinetic and boundary layer effects in the nozzle. The 
XNR2000 reactor core consists of 90 fuel elements in the outer core and 61 in the inner 
core, the equivalent inner core diameter is 29.2 cm while that of the outer core is 45.94 
cm. 

This study considers a turbulent flow of hydrogen in a circular tube with a variety of 
thermal boundary conditions. The computational models are used to evaluate and 
compare the applicability of a number of widely used correlations for Nusselt number to 
the high temperature and high heat flux conditions of hydrogen cooled nuclear rocket 



99 


2 

cores. In these systems the surface heat flux ranges from 10 to 1000 MW/m , the flow 
channel outlet temperature can be as high as 3200° K, the ratio between diameter and total 
length varies from 150 to 1000 and the Reynolds number of the flow is in the range of 
10 4 to 10 6 . The total flow length in this design is 0.55 m, and the tube diameter is 
0.0032m resulting in L/D=165. The grid mesh is 90x50. A stagnation inlet pressure of 
6.85 MPa and an back pressure of 6.55 MPa are used. The thermal boundary conditions 
and Reynolds numbers for two cases considered in this section are listed in following 
table: 

Table 5.1 Thermal conditions and Reynolds number for XNR2000 rocket core. 


Case 

Tjn (K) 

T W (K) 

WMiM MOIM 

Reynolds Number 

1 ( Uniform T w ) 

500 

500+2300z/L 

- 

20000 

2 (Uniform q”) 

500 

- 

7+1 5sin(7cz/L) 

15000 


In computation procedure, approximately 1000 time steps and 5 seconds CPU time 
per time step are required to reach a steady state. The final temperature residual is less 
than 0.00001. 

5.3.1 Non-dimensional Velocity and Temperature Profiles 

The non-dimensional velocity profiles for the uniform heat flux boundary 
condition case are presented in Figure 5-16. It can be seen that the profiles rapidly 
accelerate very near the wall before flattening out and increase slightly toward the 
centerline, which is characteristic of turbulent velocity profiles. The flow reached a fully 
developed state within 20.9 diameters of the entrance. This is very close to the results 
obtained by Barbin and Jones. [ 61 The non-dimensional temperature profiles are shown in 
Figure 5-17. This figure illustrates high temperature gradients at the inlet region and 
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flattens out downstream as the bulk temperature increases. The temperature profiles 
achieve fully developed conditions at the approximate 121.9 diameters from the tube 
entrance. These results indicate that a velocity profile achieves the fully developed regime 
before a fully developed temperature profile is obtained, that is the thermal entrance 
length did not exceed the hydrodynamic entrance length. Due to temperature dependent 
properties used for hydrogen, the temperature profiles shown in Figure 5-17 are slightly 
different for those predicted for constant properties. Viscosity and thermal conductivity of 
hydrogen are strong functions of temperature. They lead the values for hydrogen to be 
significantly different at the wall and at the centerline temperatures. 

5.3.2 Nusselt number calculations 

The Nusselt number was calculated for both cases listed in the Table 5. 1 , as shown in 
Figure 5-19. This figure indicates that the magnitude of Nusselt number appeared to be 
mainly affected by Reynolds number. The higher the Reynolds number, the higher the 
Nusselt number. The shapes of the two cases have different characteristics. Downstream 
of the entrance the Nusselt number gradually increased, after initially decreasing for 
uniform heat flux case, exhibiting a slight difference with the uniform wall temperature 
case. It is due to that the trend of Nusselt number for case 1 is mainly dependent upon the 
magnitude of dT/dr, but for case 2 it is mainly dependent upon the difference between T w 
and T b . So, these different trends between the uniform wall boundary and uniform heat 
flux boundary result in a slightly different heat transfer mechanism. 

Among the long list of correlations for the Nusselt number, four correlations which 
seem to have more relevance to the heat transfer in a nuclear rocket core are chosen for 
this analysis. It should be noted that some of these correlations have been developed 
based on experimental data which are not fully compatible with the boundary conditions 
used in cases analyzed here. These correlations are used as a point of reference not to 
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justify the validity of computational results. Figure 5-18 compares the numerical results 
with four Nusselt number correlations without axial distance and property corrections. 
These correlations include the Colburn equation, the Gnielinski correlation, the Karman- 
Boelter-Martinelli equation and the Notter-Sleicher correlation, the Gnielinski 
correlation 1161 is an updated version of the Petukhov-Kirillov correlation 1331 where the 
range of Reynolds number and Prandtl number has been expanded. For the high 
temperature and high heat flux hydrogen flow cases, a correction factor for the axial 
distance and real gas property must be used. Figure 5-19 shows a comparison between 
the numerical results and the Nusselt number correlations when the axial distance and 
property corrections are applied. The Petukhov-Kirillov-Gnielinski equation, when 
combined with the Perkins and Worsoe-Schmidt axial distance correction and the Notter- 
Sleicher property correction has the best agreement with the numerically calculated 
Nusselt number, differing by less than 3% at the last section measured. 

Figures 5-18 and 5-19 vividly show the importance of using distance and property 


corrections. 
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Figure 5-20 Comparison of the numerically calculated Nussclt number with the 
correlations without correction. 
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Figure 5-21 Comparison of the numerically calculated Nusselt number with the 
correlations including axial distance and property correction factors. 






CHAPTER 6 

SUMMARY AND CONCLUSIONS 


A computational model based on the axisymmetric, thin-layer Navier-Stokes 
equations has been proposed to investigate the convective, radiative and conductive heat 
transfer in nuclear reactors. An implicit-explicit, finite volume, MacCormack method in 
conjunction with the Gauss-Seidel line iteration procedure is utilized to solve the 
governing equations. The subsonic and supersonic flows of Hydrogen, Helium and 
Uranium Tetrafluoride under variable boundary conditions, such as adiabatic, isothermal 
and constant heat flux, are employed to simulate coolant flow on reactor cores. An 
enthalpy-rebalancing scheme is implemented to allow the convergence solutions to be 
obtained with the application of a wall heat flux. A two-dimensional method based on 
finite element technique is used to investigate the geometric behavior of nuclear reactor 
fuel elements. Using the developed models the following studies have been completed. 

1 Thin layer Navier-Stokes equations Due to the flow field at a high Reynolds 
number, the thin-layer Navier-Stokes equations are used as the governing equations in 
this study. In the thin-layer approximation to the full Navier-Stokes equations, the 
viscosity terms containing derivatives parallel to the body surface, which are in inverse of 
Reynolds number, are neglected. As a result, a substantial fraction of the available 
computer storage and time is expended in resolving the normal gradients in the boundary 
layer. In these calculations, the values of y + at nearest the wall are less than 2. So a highly 
stretched grid system is generated to lead stability and convergence. Good agreement 
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between the computed results and prior calculations in the literature indicate the validity 
of the governing equations. 

2. Numerical modeling A hybrid implicit-explicit, MacCormack scheme based on the 
finite volume approach is employed to solve the governing equations. This numerical 
algorithm has a rapid convergence compared to a simple explicit method. Normally, less 
than 1000 iterations are needed to reach the steady-state for any flow conditions and any 
fine grid system taken. MacCormack method is one of the most efficient of the second- 
order scheme and the most popular two step Lax-Wendroff method for solving problems 
with shock-capturing schemes. The calculations of this study indicate the good stability, 
consistency and efficiency. 

T Turbulence model The Baldwin and Lomax two-layer algebraic turbulence model 
is used in this study. For this turbulence model, the effects of turbulence are simulated in 
terms of an eddy viscosity coefficient p t . Thus, in stress terms of Navier-Stokes 
equations, the molecular coefficient of viscosity p is replaced by p+p t . This model 
simplifies the calculations of turbulent kinetic energy and eddy diffusivity of energy and 
avoids the necessity for finding the edge of the boundary layer. It is used to yield faster 
convergence at reasonable accuracy. 

4 Radiative beat transfer model The Rosseland diffusion approximation is used to 
simulate compressible, turbulent and developing flow of a radiative nongray gas in gas 
core reactors. A well-known Rosseland mean absorption coefficient is adopted for the 
approximation, which successfully made use of the diffusion theory in the radiation 
calculation. This method is simpler to formulate and calculate in radiation heat transfer. 
Good agreement between the computed results and the previous calculation in the 
literature indicates the accuracy of the Rosseland radiation model. 
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5. Thermal conduction model A Finite Element code ANSYS is used to investigate 
the geometric behavior of fuel elements of a reactor core, which include plate, pin, 
tubular and square lattice honeycomb. Based on the investigation a dimentionless 
number, 6 number, is defined with respect to void parameter a. The 5 number can be 
used to control temperature gradients, which is important for the lifetime of reactor 
components. The 5 number also provides an efficient method to calculate the centerline 
temperature, instead of the time-consuming CFD procedure. 

6. Wall heat transfer model and enthalpy rebalancing scheme Two boundary 

conditions, constant wall temperature and constant wall heat flux, on the outside surface 
of straight tube were included within the study. The relation of the boundary conditions is 
found based on the Fourier’s law. The difficulties in obtaining converged solutions with 
heat flux boundary conditions were addressed. An enthalpy rebalancing scheme was 
developed. The implementation of the enthalpy balancing scheme allowed convergent 
solutions to be obtained with the application of a wall heat flux. This method has been 
proven successful when problems with heat flux boundary condition are to be solved. 

7. Validation of modeling The numerical models were used to calculate the 
separation angle and drag force over a sphere which is bounded in a cylindrical tube. For 
this classic project, good agreement between numerical results and the measured data in 
the experiment indicated the validity of the developed models. Many different Nusselt 
number equation, property corrections and axial distance corrections were investigated to 
qualify the numerical calculation. Results demonstrated that the Petukhov-Kirillov- 
Gnielinski equation combined with Perkins and Worsoe-Schmidt axial distance 
correction and the Notter-Sleicher property correction is the most relevant to this study. 
This results also implied the assessment of the numerical models. 
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8. Internal heat generation The effect of internal heat generation on the heat transfer 
in the gas core reactors is examined for a variety of power densities, lOOW/cc, 500W/cc 
and lOOOW/cc. The maximum temperature, corresponding with the heat generation rates, 
are 2150°K , 2750°K and 3550°K, respectively. This analysis shows that the maximum 
temperature is strongly dependent on the value of heat generation rate and also indicated 
that a heat generation rate higher than lOOOW/cc is necessary to maintain the gas 
temperature at about 3500°K, which is typical temperature required to achieve high 
efficiency in the gas core reactors. 

9. Convective and radiative heat flux The convective and radiative heat fluxes are 
predicted for the gas core reactors. The maximum value of heat flux occurs at the exit of 
the reactor core. Radiative heat flux increases with higher wall temperature. This 
behavior is due to the fact that the radiative heat flux is strongly dependent on wall 
temperature. This study also found that at temperature close to 3500°K the radiative heat 
flux is comparable with the convective heat flux. 

The above accomplishments clearly demonstrate that the numerical model for 
convective, conductive and radiative heat transfer has the ability to simulate the heat 
transfer performance of nuclear reactor cores. This model can be used to predict the 
important information about heat transfer in nuclear reactor systems, and to qualify the 
classic empirical correlations which have been used for many years. There is no doubt 
that the numerical model is a reliable computational tool for heat transfer analysis in 
nuclear reactor analysis. 
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